In [9]:
import numpy
import matplotlib
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
%matplotlib inline
from pylab import *
In this example we want to solve the 2D example of module 4 of our class as the last module which we learned in the class. We solved this 2D problem with explicit and implicit methods. Now we want to compare these two methods (LBM and FDM)
In [10]:
nx=40 # the number of nodes in x direction lattice direction
ny=40 # the number of nodes in y direction lattice direction
alpha=0.25 # heat diffusion coefficient
D=3.0 # the dimension of the problem
mstep=100 # the number of time step
omega=1./(0.5+(alpha*D)) #Chapman-Enskog dt=1 and dx=1
Tleft=100.0 # left wall temperature
Tbottom=100.0 # right wall temperature
k=9 # k=0,1,2,3,4,5,6,8,9
In [11]:
x=numpy.linspace(0,1,nx+1)
y=numpy.linspace(0,1,ny+1)
w=numpy.ones(k) # witghting factor
T=numpy.ones((ny+1,nx+1) ) # Temperature matrix
f= numpy.ones((k, ny+1,nx+1)) # distribution function
In [12]:
w[0]=4./9 #4./9
w[1:5]=1./9.
w[5:9]=1./36.
In [13]:
T[0:ny+1,0:nx+1]=0.0
T[0:ny+1,0]=100.0
T[0,0:nx+1]=100.0
for i in range(nx+1):
for j in range(ny+1):
for l in range (k): #k=0,1,2,3,4
f[l,j,i]=w[l]*T[j,i]
In [14]:
## Main loop : comprised two parts :collision and streaming
##=====================
for n in range(mstep) :
# collision process
for i in range(nx+1):
for j in range(ny+1):
for l in range (k):
feq=w[l]*T[j,i]
f[l,j,i]=(1.-omega)*f[l,j,i]+omega*feq
#streaming process
# ==========================
for i in range(nx):
for j in range(ny,0,-1):
f[6,j,i]=f[6,j-1,i+1]
f[2,j,i]=f[2,j-1,i]
for i in range(nx,0,-1):
for j in range(ny,0,-1):
f[5,j,i]=f[5,j-1,i-1]
f[1,j,i]=f[1,j,i-1]
for i in range(nx,0,-1):
for j in range(0,ny):
f[8,j,i]=f[8,j+1,i-1]
f[4,j,i]=f[4,j+1,i]
for i in range(0,nx):
for j in range(0,ny):
f[7,j,i]=f[7,j+1,i+1]
f[3,j,i]=f[3,j,i+1]
# Boundary condition
# =============================
for j in range(0,ny+1) : #left Boundary
f[1,j,0]=( Tleft*(w[1]+w[3]) )-f[3,j,0]
f[5,j,0]=( Tleft*(w[5]+w[7]) )-f[7,j,0]
f[8,j,0]=( Tleft*(w[8]+w[6]) )-f[6,j,0]
for i in range(0,nx+1): #bottom Boundary
f[2,0,i]=( Tbottom*(w[2]+w[4]) )-f[4,0,i]
f[5,0,i]=( Tbottom*(w[5]+w[7]) )-f[7,0,i]
f[6,0,i]=( Tbottom*(w[6]+w[8]) )-f[8,0,i]
for i in range(1,nx): # top Boundary
for l in range(k):
f[l,ny,i]=f[l,ny-1,i]
for j in range(1,ny): # right Boundary
for l in range(k):
f[l,j,nx]=f[l,j,nx-1]
#================================
for i in range(nx+1):
for j in range(ny+1):
sum=0.0
for l in range (k):
sum=sum+f[l,j,i]
T[j,i]=sum
T[0:ny+1,0]=Tleft
T[0,1:nx]=Tbottom
T[1:ny,nx]=T[1:ny,nx-1]
T[ny,1:nx]=T[ny-1,1:nx]
#==============================
#==============================
In [15]:
# Finite Difference Method
#=============================
Tf=numpy.ones((ny+1,nx+1) ) # finite difference Temperature matrix
To=numpy.ones((ny+1,nx+1) )
lx=0.01
ly=0.01
x=numpy.linspace(0,lx,nx+1)
y=numpy.linspace(0,ly,ny+1)
mstep=100
dx=lx/nx
dy=ly/ny
alpha=0.25
dt=0.2*(dx**2)/alpha
print dt
Tf[0:ny+1,0:nx+1]=0.0
Tf[0:ny+1,0]=100.0
Tf[0,0:nx+1]=100.0
for n in range (mstep):
To[0:ny+1,0:nx+1]=Tf[0:ny+1,0:nx+1]
for j in range (1,ny):
for i in range (1,nx):
Txx=(To[j,i+1]+To[j,i-1])
Tyy=(To[j+1,i]+To[j-1,i])
Tf[j,i]=To[j,i]+ (dt*alpha*(Txx+Tyy-4.*To[j,i])/(dx**2))
#
Tf[1:ny,nx]=Tf[1:ny,nx-1]
Tf[ny,1:nx]=Tf[ny-1,1:nx]
##
In [16]:
title('Lattice Boltzmann Method')
CS1 = plt.contourf(x,y,T,20)
plt.colorbar();
In [17]:
title('Finite Difference Method')
CS2 = plt.contourf(x,y,Tf,20)
plt.colorbar();